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Abstract 

In this paper we describe the numerical analysis un- 
derlying our efforts to develop an accurate and re- 
liable code for simulating flame propagation using 
complex physical and chemical models. We discuss 
our spatial and temporal discretization schemes, 
which in our current implementations range in order 
from two to six. In space we use staggered meshes 
to define discrete divergence and gradient operators, 
allowing us to approximate complex diffusion opera- 
tors while maintaining ellipticity. Our temporal dis- 
cretization is based on the use of preconditioning to 
produce a highly efficient linearly implicit method 
with good stability properties. High order for time 
accurate simulations is obtained through the use of 
extrapolation or deferred correction procedures. We 
also discuss our techniques for computing stationary 
flames. The primary issue here is the automatic gen- 
eration of initial approximations for the application 
of Newton's method. We use a novel time-stepping 
procedure, which allows the dynamic updating of 
the flame speed and forces the flame front towards 
a specified location. Numerical experiments are pre- 
sented, primarily for the stationary flame problem. 
These illustrate the reliability of our techniques, and 
the dependence of the results on various code param- 
eters. 


1 Introduction 

Comprehensive mathematical models of reacting 
gases have been known for some time. (See, e.g., 
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[18].) Nonetheless, due to the complexity of the 
models, it is necessary to use numerical methods to 
extract detailed predictions and evaluate their abil- 
ity to faithfully reproduce experimental results. In 
this note we review our ongoing efforts to design and 
implement reliable numerical methods for the solu- 
tion of these complex models. Of course, we are not 
the first researchers to attempt to solve these equa- 
tions. Of particular note is the work of groups at the 
Ballistic Research Laboratory [5], Sandia National 
Laboratory [11] and the Naval Research Laboratory 
[10]. Nonetheless, some of our approaches differ from 
those pursued by other groups and seem worthy of 
additional investigation. In the design process, we 
have been motivated by potential efficiency in mul- 
tidimensional simulations. In this work, however, 
we will restrict ourselves to examples in one space 
dimension. 


2 Models 


Our starting point is the equations for conservation 
of mass, momentum, energy and species: 

^ + V • (pu) = 0, (2.0.1) 

? + (u • V)u + -V<£ = -V ■ J v , (2.0.2) 
at p p 

rs IP 

^ + (u * V)e = — V • q V • u, (2.0.3) 

at p p 

ftY 1 

V + (tx • V)Yi = — V • (pYiVi) + . (2.0.4) 

at p 

Here, p is the density of the gas mixture, u is the 
velocity vector, e is the internal energy and Yi, i = 
1, . . . , n 5 , is the mass fraction of the ith species. 

We consider a zero Mach number model, akin to 
the incompressible equations for fluid flow. Pre- 
cisely, our gas law equation of state takes the form: 


e= % ej W' (m5> 


t= 1 


i= 1 
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ei(T) = hi(T)-RT. (2.0.6) 

That is, we assume that the pressure, to leading or- 
der, is spatially const ant. The quantity <j> is the 
spatially varying pressure perturbation. It is only 
retained in the momentum equation. In a closed do- 
main, an additional equation is needed to specify the 
temporal evolution of P. In an open domain, how- 
ever, we may assume that P is constant. Detailed 
mathematical analyses of these equations and their 
derivation are given in [13, 8]. 

There are both advantages and disadvantages to 
using the zero Mach number system. The primary 
advantage is in the design of time stepping strate- 
gies, as large imaginary eigenvalues of the spatial 
operators are suppressed. Disadvantages are its re- 
stricted applicability to problems where the Mach 
number remains small and the necessity to solve ad- 
ditional equations to simulate sound production. 

To complete the problem specification, we must 
describe the computation of the reaction rates, r*, 
the diffusion velocities, V* , the heat flux, q , the en- 
thalpies, hi, and the viscosity coefficients. Our pro- 
grams are designed to accommodate any reasonable 
specifications of these quantities without requiring 
basic changes in the algorithm. In this work the 
species production rates were computed using the 
procedures (and rate coefficient expressions) built 
into the NASA Lewis kinetics code LSENS [16]. The 
NASA Lewis polynomial expressions [9] were used 
to calculate the thermodynamic properties, and the 
routines for the transport properties were adapted 
from the Sandia 1-D flame code [11]. The species 
diffusion velocities were obtained by solving the mul- 
ticomponent diffusion equation, which is a matrix 
equation — see equation E18 in [18]. 

3 Plane Flames 

A fundamental problem of computational combus- 
tion is to determine flame speeds and species profiles, 
or to establish the nonexistence of the flame, given 
the mixture composition, temperature and pressure. 
Mathematically, the plane flame is a traveling wave 
solution of the governing equations. That is, if x is 
the coordinate normal to the flame front and v is 
the flame speed we seek a solution of (2.0.1)-(2.0.4) 
which is a function of z = x + vt. Equivalently, 
we may make a Galilean transformation and seek a 
steady solution depending on x alone with a nonva- 
nishing x-velocity, u — u, at x = — oo. 

Generally, boundary conditions for traveling wave 
problems involve specifying limiting rest states at 
±oo. This description does not strictly apply to the 


plane flame problem, as there is no rest state as- 
sociated with the unburnt gas. (This is sometimes 
called the cold boundary difficulty.) Here we sim- 
ply impose an unburnt state, that is fix the unburnt 
fuel mixture, Y^o, temperature, T 0 , and pressure, P 0 . 
This raises doubts concerning the strict existence of 
the unbounded domain solution. However, so long 
as the unburnt state is an approximate equilibrium 
(that is the reaction rates are sufficiently small), we 
expect that our truncated problems do possess solu- 
tions which correctly model the physics. 

True equilibria, which will define the burnt state, 
are defined by the n s equations: 

r i (Y 1 , b ,...,Y n ^T b )= 0. (3.0.7) 

As the number of unknowns exceeds the number 
of equations by one, it is reasonable to expect 
that a one-parameter family of equilibria exist. To 
uniquely determine the burnt state, we note that un- 
der the zero Mach number assumption, an additional 
conservation law for the enthalpy can be derived: 

q + ^ p(u *f v + V^Yih^j = 0. (3.0.8) 

(Here h{ is the enthalpy of species i.) Integrating 
this equation from the unburnt to the burnt state, 
imposing the vanishing of the gradients and, hence, 
of q and V{ in these states, and taking into account 
mass conservation, we arrive at the additional rela- 
tion: 

n a ns 

Y Yi )0 ht,o = Y iyb h iyb . (3.0.9) 

i = 1 i=l 

Supplementing (3.0.7) by (3.0.9), the burnt state is 
finally determined. In practice, the chemical equilib- 
rium state is not computed using equations (3.0.7) 
and (3.0.9). It is instead obtained by minimizing the 
Gibbs function — for details see [9, 16]. 

An additional property of traveling wave solutions 
is their translation invariance. That is, given any 
solution a new solution can be obtained by adding 
any constant to the independent variable, effectively 
translating the profile. Mathematically, this implies 
that the space derivative of the solution is a nullfunc- 
tion of the linearized equations. To fix our profile 
in space we impose an additional phase condition. 
Many such conditions are possible. (For an interest- 
ing discussion of the discretization theory of travel- 
ing wave problems see Beyn [4].) What is essential is 
that the additional condition be nonvanishing when 
applied to the derivative of the solution. We choose 
a temperature normalization: 

T{x h ) = T h , To < T h < T b . (3.0.10) 
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For this to be effective we must have: 

^(x h ) > 0. (3.0.11) 

So long a s T = Th occurs within the flame front, 
(3.0.11) is expected to be valid, and has been in all 
our experiments. Later we will discuss a new method 
for dealing with the phase condition in the context of 
a time stepping procedure for computing the flame. 


4 Spatial Discretization and 
Adaptivity 

We approximate spatial derivatives using second to 
sixth order central differences, switching to lopsided 
schemes where necessary near the boundaries. We 
expect that higher order methods will prove to be 
more efficient, particularly when the physical mod- 
els are complex. In these cases, the overhead associ- 
ated with the use of higher order differences is small 
compared with the cost of evaluating various phys- 
ical quantities at each point. Hence, any reduction 
in the number of mesh points results in a reduction 
of computing time. 

Two distinct spatial difference formulas are em- 
ployed. The first, which is denoted D g and used to 
approximate first order terms in the equation, takes 
a function defined on a uniform grid and produces 
approximate derivatives on that grid. The second 
operator, Dh, approximates derivatives midway be- 
tween gridpoints. These operators are composed to 
approximate second order derivatives on the grid. 
For example, we approximate the species conserva- 
tion equations (2.0.4) by: 

uD 3 Yi = -~D h ((hpXhYJVi) + r i} (4.0.12) 
P 

where Ih is an interpolation operator of appropri- 
ate order producing function values midway between 
gridpoints from values on the grid. The diffusion ve- 
locities are computed, again midway between grid- 
points, by: 

V { = Vi {D h Yj , D h T> I h Yj, I h T). (4.0.13) 

The reason for staggering the grids in this way is 
to guarantee the ellipticity of the resulting discrete 
diffusion operators. The only grid functions anni- 
hilated by Dh are constant, whereas D g also anni- 
hilates the function (-l) j where j is a grid index. 
Therefore, the Fourier transforms of second order 
operators formed by composing Dh with itself are 
always large for large wavenumbers while those ob- 
tained by composing D g with itself decrease to zero 


as the wavenumber approaches 7r/Arc, that is as we 
approach the largest wavenumber representable on 
the grid. 

For the one-dimensional problems considered 
here, the essential length scales are the computa- 
tional domain length, L, which we cannot make too 
small without using far more sophisticated bound- 
ary treatments, and the width of the flame zone. 
In practice the ratio of these is sufficiently small 
to merit an adaptive mesh strategy. Currently, we 
use a parametrized family of coordinate maps. This 
strategy has been successfully employed and ana- 
lyzed by Bayliss and coworkers [1, 2] in the context 
of spectral simulations of simpler combustion mod- 
els. The maps we use are slightly different. Precisely, 
we write: 


1 

x = xi + -7 y 



(1 -e)(y/0y \ 

i + (y/P) p J ’ 


(4.0.14) 


and employ a uniform grid in y. For e < 1 there is a 
clustering of gridpoints near x = xi. Using the tem- 
perature as our monitor function, we choose xi to 
coincide with the location of the maximum deriva- 
tive of T and choose c according to the ratio: 


Tt-Tp 

e L(T x ) max 


(4.0.15) 


Throughout our simulations we have chosen p = 8. 
The parameters 7 and (3 are chosen to scale the 
length of the y-interval and to guarantee that suffi- 
ciently many points are in the refinement layer. Our 
exact formulas are: 

/3 = (l + 3e)- 1 , 7 = ^^ L. (4.0.16) 

A defect of the very simple mapping described 
above is its inability to handle multiple fronts. An 
alternative method we have used which is capable of 
handling multiple fronts is given by the mapping: 

y = aJ* (1 + a\T x \ + (1 - 01 T xx \ 1/2 ) 1/2 dx. 

0 (4.0.17) 

where A is a normalization parameter and 0 < f < 1. 

In all cases we use static rezoning, that is update 
the mesh after some fixed number of time steps and 
interpolate the solution onto the new mesh. Each 
of these methods would require extensive modifica- 
tion for use in multiple dimensions. Therefore, a 
more flexible approach such as the AMR techniques 
of Berger and coworkers (e.g. [3]) will be considered. 


5 Temporal Discretization 

It is well known that the equations of reacting flows 
may be stiff when realistic reaction mechanisms are 
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included. Therefore, if an explicit temporal dis- 
cretization were employed, the time step could be 
severely restricted by stability requirements. Fully 
implicit discretizations, on the other hand, require 
the solution of a discrete approximation to a nonlin- 
ear elliptic system at each time step. Solving such 
a problem is potentially quite expensive in multiple 
dimensions. Here we take an intermediate approach, 
based on the concept of preconditioning. 

Consider, for simplicity, a system of ordinary dif- 
ferential equations: 

- = F{u). (5.0.18) 

Introduce a preconditioning matrix, G(u). A first 
order approximation to (5.0.18) is given by: 

(7 - A tG(u(t)))w = AtF(u(t)), (5.0.19) 

u(t + At) — u(t) + w . (5.0.20) 

For reasons of efficiency, we want to choose G so that 
(5.0.19) may be easily solved. To study the stability 
of the approximation, we consider the special case 
F = Ju. Then we have: 

u(t+At) = (/-AtC?) -1 (I+At(J-G))u(t) (5.0.21) 


order time advancement procedures based on the 
preconditioning strategy: extrapolation, as in the 
ode solver SEULIM (e.g. [6]), and iterated deferred 
correction, as discussed in [7]. As the second ap- 
proach is less familiar, and is the one employed in 
the example given below, we describe it briefly. 

Suppose u(t n ) is our approximate solution at t = 
t n and that we wish to compute u(t n + 1 ). Recall the 
formula obtained by integrating (5.0.18) from t n to 
an arbitrary time, t: 

u(t) = u(t n ) + [ F(u($))ds. (5.0.25) 

Jt n 

Choose nodes (typically of Gauss, Radau or Lobatto 
type), 

t n ^ Si <1 S 2 ^ . . . <C $m ^ tn-j-i- (5.0.26) 

Let the pth approximate solution on the nodes be 
given by: 

v] « u( Sj ), (5.0.27) 

where the inital approximation is computed by 
(5.0.19), (5.0.20) with At replaced by — Sj. To 
compute corrections, define the interpolant, Q p (s), 
to (sj,F(Vj)) and the residual: 


so that the method is stable if and only if for some 
subordinate matrix norm: 



Q p (s)ds - v p . 


||(7 - AtG)~ l (I + At(J - G))|| < 1. (5.0.22) Then solve: 


(5.0.28) 


As a simple example, suppose J is a discrete approxi- 
mation to a variable coefficient multidimensional dif- 
fusion operator: 

J = '^2D h ,j(cr j (x)D h j) i (5.0.23) 

j 

where Dh are as discussed above. Then it is possible 
to take: 

(7 - AtG) = n<' -A t&jD+jD-j) (5.0.24) 

j 

for suitably chosen oj. Then we only need to solve 
small tridiagonal systems, no matter what stencil 
width our actual difference approximations employ. 

The time stepping scheme, as outlined above, is 
only first order. This is sufficient for the compu- 
tation of stationary flames, where time-stepping is 
simply used to generate an initial approximation for 
Newton’s method. For time accurate simulations, on 
the other hand, it is more efficient to use a higher 
order method. There are two distinct approaches we 
have considered for automatically producing higher 


(I ~ kjG(vj))(Sj +1 — Sj) =(5.0.29) 

hj — Sj , (5.0.30) 

t^ +1 = V? + f%. (5.0.31) 

The overall method order increases by one at each 
iteration, up to a maximum determined by the inte- 
gration scheme. Assuming this order is not less than 
Z, an 1th order approximation to u(t n+ i) is given by: 

u(t n+ 1 ) = u(* n ) + f n+I Q t - 1 (s)ds. (5.0.32) 

Note that our semidiscretized system is in fact a 
differential-algebraic rather than a differential equa- 
tion. Formally, the index of this system is three. 
In particular, no equation for the the time deriva- 
tive of the pressure variation, 0, is given. To obtain 
such an equation, one must differentiate the con- 
straint, namely the equation of state (2.0.5), three 
times. (The first differentiation produces pt] substi- 
tuting for p t using (2.0.1), the second differentiation 
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produces u t \ finally, substituting for u t using (2.0.2) 
yields an equation for <f> which can be differentiated 
and solved for <j> t .) In one space dimension, however, 
(2.0.2) and, hence, <f> may be eliminated. Then we 
have no equation for ut, and the index, following the 
same reasoning as above, is reduced to two. 

In the current code, we only apply our time- 
stepping procedure to the spatially discretized ap- 
proximations to (2.0.3) and (2.0.4). The density is 
computed using (2.0.5) and the velocity is updated 
using (2.0.1). For stationary flames this reduces to: 

u-^-v, (5.0.33) 

P 

where v is the velocity at x = — oc. For time ac- 
curate approximations we found it difficult to use 
(2.0.1) directly. Instead, we replace p t by first ex- 
pressing it in terms of e t and Y^t and then using 
(2.0.3) and (2.0.4) to eliminate all time derivatives. 
This results in a somewhat complicated expression, 
which is nevertheless easily integrated for u. 

In our code, G is chosen by replacing the heat 
conduction and species diffusion terms in equations 
(2.0.3) and (2.0.4) by simplified physical approxima- 
tions and by using second order differences. It also 
includes the exact Jacobian of the reaction rates. 

Finally, we note that the time steps are auto- 
matically adjusted, within preset upper and lower 
bounds, according to the relative change in the so- 
lution. 

6 Stationary Flames 

As discussed in Section 3, the plane flame prob- 
lem is to compute stationary solutions to equations 
(2.0.1)-(2.0.4) along with the normalization condi- 
tion (3.0.10) and the boundary conditions: 

(e, Yi) -+ (eo,Fi,o), * -> -oo, (6.0.34) 

(e,Yi) -4 {e b ,Y itb ), x -4 oo, (6.0.35) 

u — > v, x —oo. (6.0.36) 

Recall that the state (e 0 ,li,o) is specified, the state 
(e b ,Yi 9 b) is computed given (e 0 ,i<,o) by minimizing 
the Gibbs function, and the flame speed, t/, is to be 
detrmined as part of the calculation. 

A direct approach to the computation is to trun- 
cate the infinite domain by a finite domain, replac- 
ing (6.0.34)-(6.0.36) by approximate boundary con- 
ditions at the artificial boundaries, then to approxi- 
mate the spatial derivatives by differences on a (gen- 
erally nonuniform) mesh, and finally to use some 
variant of Newton’s method to solve the resulting 


algebraic system. The problem here is in the last 
step, as the convergence of Newton’s method is only 
guaranteed for a sufficiently accurate initial guess, 
which seems difficult to automatically generate for 
this problem. A more reliable alternative is to use 
a time-stepping procedure in place of Newton iter- 
ations. In this approach one makes use of the ap- 
parently large dynamical basin of attraction of the 
plane flame solution. Disadvantages associated with 
standard time-stepping procedures include slow con- 
vergence to steady state, an inability to fix the flame 
front, and difficulties in determining a precise speed. 

In this work we use a hybrid approach. We begin 
with a novel time-stepping procedure, which allows 
us to dynamically compute a flame speed and to 
modify the dynamics to push the front towards its 
normalized location. Similar ideas have been pro- 
posed by Sermange [17], but ours differ in the details. 
Second, either after the time derivatives become suf- 
ficiently small or after a maximum number of time 
steps are carried out, we switch to a damped Newton 
iterator to achieve final convergence. 

The time-stepping procedure is most easily ex- 
plained in a general setting. Consider the nonlinear, 
algebraic eigenvalue problem: 

F(w) - vSw = 0. (6.0.37) 

Here v is the eigenvalue and 5 is a square matrix. 
Discrete approximations to the plane flame problem 
take this form when the velocity, u, is replaced using 
(5.0.33) and p is eliminated using (2.0.5). Equation 
(6.0.37) then represents the discrete approximation 
to (2.0.3)-(2.0.4) multiplied by p. The matrix S is 
the discrete derivative operator. 

As mentioned above, solutions to the plane flame 
problem are only unique up to translation. This 
implies that the equations linearized about a plane 
flame solution have a nullfunction which is the 
derivative of the solution itself. Moreover, we de- 
termine a unique solution by imposing an additional 
normalization condition (3.0.10). We represent our 
discrete approximation to (3.0.10) by 

N(w) = 0, (6.0.38) 

where N is scalar valued. Corresponding to the 
translation invariance, we assume that if w* is the 
solution of (6.0.37)-(6.0.38), then the matrix: 

A = DF(w*) - vS, (6.0.39) 

has 0 as a simple eigenvalue with normalized eigen- 
vector q 0 = StuV||Su/*||. (This assumption is gen- 
erally false after domain truncation. However, drop- 
ping it complicates the analysis beyond the scope of 
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this note.) We further assume that is dynami- 
cally stable, that is that all other eigenvalues, A j, of 
A, have negative real part. Moreover, we impose the 
nondegeneracy condition: 


Finally, we consider the boundary conditions im- 
posed at the artificial boundaries introduced to trun- 
cate the infinite domain. At the cold boundary we 
simply impose the Dirichlet conditions: 


DN{w*)Sw* > 0. (6.0.40) 


(e,li) — (eo,l{,o)- (6.0.49) 


Note that if the left-hand side of (6.0.40) were zero, 
the normalization would fail to guarantee a unique 
solution. In our case, (6.0.40) is equivalent to 
(3.0.11). 

Our time-stepping procedure is to replace (6.0.37) 
by the dynamical problem: 


w t — G{w ), (6.0.41) 


with G defined by: 


G(w) = F(w) — ($(w) + fj,N(w))Sw , (6.0.42) 

where $ is a scalar function and fx > 0 is a damping 
parameter. In order for w* be a rest point of the 
new system, it is necessary to guarantee that: 


$(u;*) = v. 

This we accomplish by choosing: 

= (StvfRFjtv) 

{ 1 ( Sw) T R(Sw ) ’ 


(6.0.43) 

(6.0.44) 


for some symmetric semidefinite matrix R. In our 
current code, R is simply the projection onto the 
internal energy. 

The effectiveness of our approach depends on the 
spectrum of DG(w*), which must lie in the left half- 
plane. By direct computation we find: 


DG(w*) = (/-(&$$§£§) A (6-0.45) 

— fi(Sw*)DN(w *). 


Let Q be an orthogonal matrix whose first column 
is 4o- Then we have: 

Q TAQ = ( o jl) ) . (6-0.46) 

and 

Q t DG(w*)Q = ( A 0 ° ) ’ ( 6 - 0 - 47 ) 

A 0 = — / aDN(w*)Sw* < 0. (6.0.48) 

Therefore, the spectrum of DG(w*) is given by the 
eigenvalues of A^ x \ which by (6.0.46) are also the 
nonzero eigenvalues of A, and Ao. Hence, the local 
dynamic stabilty of w* is assured. 


Our justification for this is the rapid exponential de- 
cay to this state which we observe in the compu- 
tations. At the hot boundary, where the observed 
decay rates are much slower, we impose the condi- 
tions analyzed by Loheac [12] for advection-diffusion 
equations: 

e** = Vi, xx = 0. (6.0.50) 

Although our experiments indicate that these con- 
ditions do provide sufficient accuracy, more sophis- 
ticated conditions could be used and might lead to 
improvements. See [4] for a detailed discussion. 


7 Results 


In this section we present the results of some numer- 
ical experiments with our code. The vast majority 
of these involve computations of a stationary, sto- 
ichiometric hydrogen-air flame at a temperature of 
300 degrees Kelvin and atmospheric pressure. Our 
mechanism includes 13 species. 

This problem is rather standard, and was included 
in a series of benchmark problems in laminar flame 
computations sixteen years ago [15]. We note that 
our mechanism does differ from the one used there. 
Our complete set of results is summarized in Table 1. 
Our computed flame speeds vary between 235.77cms 
and 239.08cms. If the somewhat less reliable sixth 
order results are ignored, this variance is reduced 
from 236.86 to 237.92 - i.e. about one-half of one 
percent. It is interesting to note that these results 
show far better agreement with the experimental re- 
sults of Milton and Keck [14] than do any presented 
in [15], though this no doubt says more about the 
accuracy of our models than it does about the merits 
of our code. Plots of the temperature, velocity, and 
various mole fractions for a typical case (fourth or- 
der, 149 points, domain width L = 1) are presented 
in Figures 1 and 2. 

In all cases we tried except one, the code success- 
fully ran to completion. (The one failure was with 
the sixth order differencing and our coarsest mesh, 
n = 49.) The maximum number of time steps was 
set to 7000 and the code automatically switched to 
damped Newton iterations whenever: 


\\dw/dt\\ 

IMI 


< 5sec 1 . 


(7.0.51) 
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(Here w is the complete discrete solution.) The con- 
vergence criterion for Newton’s method was a rela- 
tive change less than 10“ 5 . This criterion was not 
always achieved, in which case the final estimated 
relative error is listed in the tolerance column. Such 
failures to fully converge are most likely due to a 
highly ill-conditioned Jacobian. Note that Newton’s 
method fails to converge at all if the initial con- 
ditions for the time-stepping are used as its initial 
guess. 

Mesh Refinement: Results are shown for varying 
numbers of mesh points and various orders. Fixing 
all other parameters, the results clearly converge un- 
der refinement. 

Varying Order: Both the second and fourth order 
codes show consistent convergence behavior. How- 
ever, fixing all other parameters, the flame speeds 
computed using the second order code are slightly 
smaller. One would, of course, hope for even better 
agreement here. Again this may be attributable to 
conditioning problems. The results for the sixth or- 
der code are a little worse, at least as judged by its 
convergence under mesh refinement. It is possible 
that the sixth order method is more susceptible to 
stability problems due to the large one-sided stencils 
it uses near the boundaries. It may be necessary to 
use some additional mesh clustering at the boundary 
to improve its performance. 

Varying Domain Length, L : We have computed 
solutions on both a small domain, L = 1cm, and 
a large domain, L = 10cm. Again, for the second 
and fourth order codes there is excellent agreement 
(about 1/20 of one percent) between the results on 
diff erent domains. This indicates that our artificial 
boundary conditions are adequate for the domain 
lengths considered. 

Varying x^: Fixing our domain to be x € [0,1], 
we used the fourth order method with 149 points 
to compute solutions with varying from 1/4 (our 
standard value) to 3/4. The flame speed varied be- 
tween 236.86 and 237.90. Results for Xh € [1/4, 1/2] 
showed far less variance. We expected this, as the 
decay rate of the solution to the unburnt state is far 
greater than its decay to the burnt state. 

Varying p: We also varied the damping parameter, 
//, with the fourth order method, xh = 1/4, L — 1 
and n = 149. There were significant changes in the 
number of time steps required, with the best results 
obtained with fi — 1000. Surprisingly, the flame 
speeds also varied somewhat. Note that varying 
only changes the initial guess from the point of view 
of the Newton iterator. The discrepancies in the 
results are another indication of ill-conditioning. 
Time Accurate Simulations: We have, at this 


time, far less experience running our code for time- 
dependent problems. Therefore, we include only one 
preliminary result. It involves using the stationary 
flame profile as initial data and setting the inlet ve- 
locity to zero. The flame then propagates to the 
right. From the graphs it is clear that the pro- 
file and speed are correctly maintained. The data 
shown were obtained using the fourth order method 
in space and time with 99 mesh points. 

8 Conclusions 

Our fundamental conclusion, based not only on the 
results shown here, but on numerous others to be 
published at a later date, is that our code is both ac- 
curate and robust for stationary plane flame calcula- 
tions. We must, of course, carry out many more ex- 
periments to assess its behavior for time-dependent 
and multidimensional problems. The results here do 
point to some difficulties in achieving high accuracy, 
which are most likely due to poor conditioning. This 
issue should be investigated further. 
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Appendix 


Order 

i 

Mesh Pts. 


M 

Time Steps 

Tol. 

V 

2 

1 

49 

L/4 

50 

1514 

9 x 10" 4 

237.35 

2 

1 

49 

L/4 

50 

1514 

9 x 10~ 4 

237.35 

2 

1 

74 

L/4 

50 

2342 

1 x 10- 5 

237.21 

2 

1 

99 

L/4 

50 

1022 

1 x lO"* 

237.32 

2 

1 

149 

L/4 

50 

1417 

1 x 10" 5 

237.46 

2 

1 

199 

L/4 

50 

1845 

1 x 10 -5 

237.47 

2 

1 

249 

L/4 

50 

1838 

1 x 10" 5 

237.46 

4 

1 

49 

L/4 

50 

7000 

1 x 10" 3 

237.43 

4 

1 

74 

L/4 

50 

3503 

1 x 10~ 5 

237.81 

4 

1 

99 

L/4 

50 

2032 

1 x 10" s 

237.41 

4 

1 

149 

L/4 

50 

1550 

1 x 10~ s 

237.92 

4 

1 

199 

L/4 

50 

1507 

1 x 10-* 

237.82 

4 

1 

249 

L/4 

50 

1525 

1 x lO" 5 

237.78 

6 

1 

74 

L/4 

1000 

2716 

1 x 10" s 

238.52 

6 

1 

99 

L/4 

1000 

1089 

1 x lO" 5 

236.16 

6 

1 

149 

L/4 

1000 

714 

2 x 10" 4 

235.64 

6 

1 

199 

L/4 

1000 

775 

1 x 10-* 

238.61 

4 

1 

149 

L/4 

10 

2524 

1 x 10~ 5 

236.91 

4 

1 

149 

L/4 

250 

1023 

1 x 10" 5 

237.96 

4 

1 

149 

L/4 

500 

799 

1 x 10" 5 

236.86 

4 

1 

149 

L/4 

1000 

768 

1 x 10" b 

236.86 

4 

1 

149 

L/4 

2000 

903 

1 x 10" 5 

236.88 

4 

1 

149 

i/3 

1000 

780 

1 x 10~ & 

237.10 

4 

1 

149 

i/2 

1000 

759 

1 x 10" 5 

236.98 

4 

1 

149 

2i/3 

1000 

766 

1 x lO" 5 

237.90 

4 

1 

149 

3i/4 

1000 

735 

1 x 10~ 5 

237.63 

2 

10 

999 

i/4 

50 

5434 

1 x 10" 6 

237.31 

2 

10 

1999 

i/4 

50 

5352 

1 x 10" 5 

237.48 

2 

10 

2999 

i/4 

50 

5670 

1 x 10" s 

237.50 

4 

10 

999 

i/4 

50 

7000 

1 x 10" s 

236.86 

4 

10 

1499 

i/4 

50 

7000 

1 x 10~ 5 

237.92 

4 

10 

1999 

i/4 

50 

5381 

1 x 10" 5 

237.92 

6 

10 

999 

i/4 

1000 

1517 

1 x 10" 5 

235.77 

6 

10 

999 

i/4 

50 

5030 

1 x 10~ 5 

236.09 

6 

10 

1499 

i/4 

50 

5713 

1 x 10~ 5 

238.81 

6 

10 

1999 

i/4 

1000 

1722 

1 x 10" s 

236.13 

6 

10 

1999 

i/4 

50 

5313 

1 x 10~ s 

239.08 


Table 1: Summary of Results 


9 








Temperature 



Figure 3. Unsteady Sol. n = 99, L = 1. 
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at time t=0, 2d-4, 4d-4, 6d-4, 8d-4 


Figure 4. Unsteady Sol. n = 99, L = 1. 
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Mole fraction of species H2 



at time t=0, 2d-4, 4d-4, 6(3-4, 8d-4, ld-3 


Figure 5. Unsteady Sol. n = 99, L = 1. 
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at time t=0, 2d-4, 4d-4, 6d-4, 8d-4 


Figure 6. Unsteady Sol. n = 99, L = 1. 
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Mole fraction of species H20 



at time t=0, 2d-4, 4d-4, 6d-4, 8d-4, ld-3 


Figure 7. Unsteady Sol. n = 99, L = 1. 
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Figure 8. Unsteady Sol. n = 99, L = 1. 
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